{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## the fft function\n",
    "Let's make sure we're doing the fft with numpy correctly (i.e. it should agree with MatLab), but to how many decimal places?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import sys\n",
    "sys.path.append(\"D:/StanfordYearTwo/RCWA\")\n",
    "from convolution_matrices import convmat2D as cm\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.01851004+0.j         0.80316487+1.39112235j 0.80316487-1.39112235j]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([ 0.87494659,-0.73138314,0.87494659])\n",
    "\n",
    "result = np.fft.fft(A,n=3)\n",
    "print(result)\n",
    "## appears to be agreement to a large number of decimal places"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "ans = [1.018510040000000 + 0.000000000000000i, \n",
    "  0.803164865000000 + 1.391122353034198i, \n",
    "  0.803164865000000 - 1.391122353034198i]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### in RCWA\n",
    "In RCWA, for patterned layers, we can take an FFT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATEAAAD8CAYAAAAfZJO2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAFQlJREFUeJzt3X2MXNV5x/HvL2vAgbwYbIMc26mJ\nYqWgKLzIIk6pWmJCY2gU8wdI0Ly4qSWrKm1JEykxitQoUv9I1CqkkSLUbUxxWkpICQgLoTgWL4oi\nFRMDDoE4xBua4o1djHlL2ogX7z79456FYb32np2duXPO7u8jHe3cM3fuPLP2ffacc889o4jAzKxW\nbxp0AGZms+EkZmZVcxIzs6o5iZlZ1ZzEzKxqTmJmVrW+JDFJ6yU9IWlE0pZ+vIeZ1UvSjZIOSXqs\no+7vJf1M0qOS7pC0KOdYPU9ikoaAbwCXAmcDV0s6u9fvY2ZVuwlYP6luJ/DeiHgf8HPgupwD9aMl\ndgEwEhFPRsQrwLeBDX14HzOrVET8AHhuUt33I+JI2nwAWJFzrAU9jg1gObC/Y3sUeP/xXrDktKFY\ntfKEPoRiZhMeevTlwxGxtNvXf/iDp8Szz43lvtfjwEsdVcMRMTyDt/sz4NacHfuRxDRF3VH3Nkna\nDGwGeOfyBTy4Y2UfQjGzCUPLRv57Nq9/9rkxHtzxzsz32vdSRKzp5n0kfQE4Atycs38/ktgo0JmR\nVgAHJu+UsvIwwJpzFvoGTrPCBTDOeF/fQ9JG4CPAxZF5Y3c/ktiPgNWSzgR+BVwF/Ekf3sfMWhQE\nr0Zed7IbktYDnwf+MCJ+m/u6niexiDgi6S+BHcAQcGNEPN7r9zGz9vWqJSbpFuAiYImkUeCLNFcj\nTwJ2SgJ4ICL+fLpj9aMlRkTcDdzdj2Ob2WAEwViPlu6KiKunqN7azbH6ksTMbG4aP/oa3cA5iZlZ\nlgDGnMTMrGZuiZlZtQJ4tcDl7J3EzCxLEO5OmlnFAsbKy2FOYmaWp5mxXx4nMTPLJMamvDV6sJzE\nzCxLM7DvJGZmlWrmiTmJmVnFxt0SM7NauSVmZlULxFiBX5DmJGZm2dydNLNqBeKVGBp0GEdxEjOz\nLM1kV3cnzaxiHtg3s2pFiLFwS8zMKjbulpiZ1aoZ2C8vZZQXkZkVyQP7Zla9Mc8TM7Naeca+mVVv\n3FcnzaxWzQ3gTmJmVqlAvOrbjsysVhF4squZ1Uye7Gpm9QrcEjOzynlg38yqFajIRRHLS6tmVqTm\nK9sWZJXpSLpR0iFJj3XUnSZpp6R96eepOXE5iZlZpubLc3NKhpuA9ZPqtgD3RMRq4J60PS0nMTPL\nEjQz9nPKtMeK+AHw3KTqDcC29HgbcHlOXNO+20yafWp8XdKIpEclnZ8ThJnVYQYtsSWSdneUzRmH\nPyMiDgKkn6fnxJTTEruJ/GbfpcDqVDYDN+QEYWbli9BMWmKHI2JNRxnuV1zTJrEZNvs2AN+KxgPA\nIknLehWsmQ1OM7A/lFW69PREvkg/D+W8qNsxsWM1+5YD+zv2G011R5G0eaKp+cyzY12GYWbtadbY\nzyld2g5sTI83AnfmvKjXA/tTXZaIqXaMiOGJpubSxeXdVGpmb9QM7CurTEfSLcB/Au+RNCppE/Bl\n4BJJ+4BL0va0up3s+rSkZRFxcFKzbxRY2bHfCuBAl+9hZoXp1Yz9iLj6GE9dPNNjdRvRsZp924FP\npquUa4EXJ7qdZla3iRn7vWiJ9dK0LbHU7LuI5pLpKPBFmmbed1IT8CngyrT73cBlwAjwW+BTfYjZ\nzAakyi8KmUmzLyICuGa2QZlZeSLg1fEKk5iZGUx0J53EzKximfdFtspJzMyyTEyxKI2TmJllcnfS\nzCrnNfbNrFrN1cny7q5xEjOzLKUuT+0kZmbZ3J00s2r56qSZVc9XJ82sWhHiiJOYmdXM3Ukzq5bH\nxMysek5iZlYtzxMzs+p5npiZVSsCjnhRRKvBh99xbs+OtePAnp4dywbP3UkrVi8T17GO64RWN4+J\nWXH6lbimez8ns3qFk5gNWtuJa7oYnNDq4oF9G5gSktdU3DqrR4THxGxASk1gnZzMaiDGfHXS2lRD\n8prMyaxsHhOzVtSYvCb78DvOdSIrTKn3TpbXNrRZmQsJbMJc+ixzQjTjYjmlTW6JzRFz9YR397Is\nJV6ddEtsDpirCazTfPiMpYs0sJ9Tckj6G0mPS3pM0i2SFnYTl5NY5ebTyT2fPmupetWdlLQc+Gtg\nTUS8FxgCruomJiexis3Hk3o+fuaSRCirZFoAvFnSAuBk4EA3MTmJVWo+n8zz+bMPUtPKyk5iSyTt\n7iib33is+BXwD8BTwEHgxYj4fjdxOYlVyCexfweDMh7KKsDhiFjTUYY7jyPpVGADcCbwDuAUSR/v\nJiYnscr45H2dfxft6+EUiw8B/xURz0TEq8DtwO91E9O0SUzSSkn3SdqbriRcm+pPk7RT0r7089RU\nL0lflzQi6VFJ53cTmJmVJRDj42/KKhmeAtZKOlmSgIuBvd3ElfNuR4DPRsRZwFrgGklnA1uAeyJi\nNXBP2ga4FFidymbghm4Cs6O55XE0/07aFZll2uNE7AJuAx4GfkKTi4aP+6JjmDaJRcTBiHg4Pf4N\nTbZcTtOf3ZZ22wZcnh5vAL4VjQeARZKWdROcvc4n67H5d9OSmQ3sT3+4iC9GxO9GxHsj4hMR8XI3\nYc1oTEzSKuA8YBdwRkQcTMEcBE5Puy0H9ne8bDTVmfWNE1lLetUU66Hs244kvQX4LvDpiPh1042d\netcp6o76WOmS62aAdy733U/H4xPUSlHiKhZZLTFJJ9AksJsj4vZU/fRENzH9PJTqR4GVHS9fwRST\n2CJieOLy69LFQ93Gb/YaJ/v+CmB8XFmlTTlXJwVsBfZGxFc7ntoObEyPNwJ3dtR/Ml2lXEszie1g\nD2OeV3xiWjECCOWVFuW0xC4EPgGsk7QnlcuALwOXSNoHXJK2Ae4GngRGgH8G/qL3Yc8PTmAz599Z\nf1W5FE9E/JCpx7mgmdsxef8ArpllXGZd84KKfdRygsrhEfVCuUVh5ZnRzd2tcRIzs3wFtsR876TN\nSW7J9kFAjCurtMktsQL5BLRyldeddEvM5iz/MeiDmmfsm5mVOCbmJFYYtx6sWBOTXQvj7qTNaf6j\n0FtVTnY1M3tNy1cecziJFcStBiudPCZmZtUawJXHHE5iZpap/RUqcjiJmVk+t8TMrGrjgw7gaJ5i\nYXOeL5j0SKGLIrolZmbZfHXSzOpWYBJzd9LMquaWmJllc3fSzOoV+LYjM6ucW2JmVjN3J80GwF/f\n1kNOYmZWNScxM6uVoszupOeJmVm+ceWVDJIWSbpN0s8k7ZX0gW5CckvMzLL1uCX2j8D3IuIKSScC\nJ3dzECcxM8vXoyQm6W3AHwB/ChARrwCvdHMsdycL4qtoVrR4fVxsugIskbS7o2yedLR3Ac8A/yLp\nEUnflHRKN2E5idmc5j8MPZb/5bmHI2JNRxmedKQFwPnADRFxHvB/wJZuQnISM7NsGs8rGUaB0YjY\nlbZvo0lqM+YkVhi3HGw+iIj/AfZLek+quhj4aTfH8sC+zVn+g9AHvb06+VfAzenK5JPAp7o5iJNY\ngXYc2OMlla08PZ7sGhF7gDWzPc603UlJCyU9KOnHkh6X9KVUf6akXZL2Sbo1ZVMknZS2R9Lzq2Yb\npJkVIn9gvzU5Y2IvA+si4hzgXGC9pLXAV4DrI2I18DywKe2/CXg+It4NXJ/2M2uVu5J9UmMSi8b/\nps0TUglgHc0VBYBtwOXp8Ya0TXr+YknlraRWOJ+EVhrR06uTPZN1dVLSkKQ9wCFgJ/AL4IWIOJJ2\nGQWWp8fLgf0A6fkXgcVTHHPzxES4Z54dm92nMOvgPwB9MrPJrq3JSmIRMRYR5wIrgAuAs6baLf2c\nqtV11MeKiOGJiXBLFw/lxjuv+GS04tTYnewUES8A9wNrgUWSJq5urgAOpMejwEqA9Pzbged6EazZ\ndJz4+6zGJCZpqaRF6fGbgQ8Be4H7gCvSbhuBO9Pj7Wmb9Py9EVHgKkR18ElpJSmxO5kzT2wZsE3S\nEE3S+05E3CXpp8C3Jf0d8AiwNe2/FfhXSSM0LbCr+hC32VGc8FtQYHNk2iQWEY8C501R/yTN+Njk\n+peAK3sSnQGe/GqFiPavPObwvZOVcCvj+Pz7aUmNY2JWDp+oU/PvpT21jomZFcsJrGUFjom5JVYZ\nn7Q2MLldSbfEbDoe6HcyHwThr2yzHvJJbINQ4piYk1jF5msim6+fuwjuTlqvzaeupZNXAdydtH7w\nyW2tKHQVC7fE5oiJRDYXW2VO0gVxS8z6ba6d8HPt89Su2kURrS5z5cSfK59jLnF30lpTc/fSyatQ\nA7jymMNJbI6rKZk5eVXAScwGpeSpGE5edSh1xr6T2DzSmSxKSGhOXvXReHlZzElsnhpUN9OJq2Ie\nE7MStdE6c+KaO9ydtKIdL9lMl+CcqOYJJzGrlZOUgVtiZla7ApOYZ+ybWZ7o/W1HkoYkPSLprm7D\nckvMzLL0aZ7YtTRfxv22bg/glpiZ5YvIKxkkrQD+GPjmbEJyS8zMss2gJbZE0u6O7eGIGJ60z9eA\nzwFvnU1MTmJmlmdmk10PR8SaYz0p6SPAoYh4SNJFswnLSczMsvVwrbALgY9KugxYCLxN0r9FxMdn\neiCPiZlZtl5dnYyI6yJiRUSsAq4C7u0mgYFbYmaWK8getG+Tk5iZZevHjP2IuB+4v9vXO4mZWb7y\nGmJOYmaWx4simlndIopcFDH76uTke5wknSlpl6R9km6VdGKqPyltj6TnV/UndDNrXWSWFs1kisXE\nPU4TvgJcHxGrgeeBTal+E/B8RLwbuD7tZ2ZzQIlf2ZaVxCbf4yRJwDrgtrTLNuDy9HhD2iY9f3Ha\n38xqFsB45JUW5bbEJu5xmpjGthh4ISKOpO1RYHl6vBzYD5CefzHt/waSNkvaLWn3M8+OdRm+mbWq\nxu5k5z1OndVT7BoZz71eETEcEWsiYs3SxUNZwZrZYJXYncy5OnnUPU40LbNFkhak1tYK4EDafxRY\nCYxKWgC8HXiu55GbWeuqvDp5jHucPgbcB1yRdtsI3Jkeb0/bpOfvjSjwXgUzm5ncrmRp3cnj+Dzw\nGUkjNGNeW1P9VmBxqv8MsGV2IZpZCZrJrpFV2jSjya6d9zhFxJPABVPs8xJwZQ9iM7PS9G4pnp7x\njH0zy9Z2KyuHk5iZ5RnAeFcOJzEzy1TmvZNOYmaWz91JM6tW9HSN/Z5xEjOzfG6JmVnVysthTmJm\nlk/j5fUnncTMLE/gya5mVi/R/i1FOZzEzCyfk5iZVc1JzMyq5TExM6udr06aWcXC3Ukzq1jgJGZm\nlSuvN+kkZmb5PE/MzOpWYBKbzReFmNl8EgFj43llGpJWSrpP0l5Jj0u6ttuw3BIzs3y9a4kdAT4b\nEQ9LeivwkKSdEfHTmR7ISczM8vUoiUXEQeBgevwbSXuB5YCTmJn1SQD5a+wvkbS7Y3s4Ioan2lHS\nKuA8YFc3YTmJmVmmgMieY3E4ItZMt5OktwDfBT4dEb/uJionMTPLE2QN2ueSdAJNArs5Im7v9jhO\nYmaWr0djYpIEbAX2RsRXZ3MsT7Ews3wReWV6FwKfANZJ2pPKZd2E5JaYmWXq3Q3gEfFDQL04lpOY\nmeUJwEvxmFnVCrztyEnMzDJFT69O9oqTmJnlCYj8eWKtcRIzs3z5M/Zb4yRmZvkKHBPLmicm6ZeS\nfpLmcuxOdadJ2ilpX/p5aqqXpK9LGpH0qKTz+/kBzKwlEc3VyZzSoplMdv1gRJzbcT/UFuCeiFgN\n3JO2AS4FVqeyGbihV8Ga2YD1brJrz8xmxv4GYFt6vA24vKP+W9F4AFgkadks3sfMihDE2FhWaVNu\nEgvg+5IekrQ51Z2R1gSaWBvo9FS/HNjf8drRVPcGkjZL2i1p9zPPtvuhzawLE0vx5JQW5Q7sXxgR\nBySdDuyU9LPj7DvVrQRHfaq0ttAwwJpzFpY3WmhmR6t1ikVEHEg/D0m6A7gAeFrSsog4mLqLh9Lu\no8DKjpevAA70MGYzG4AAosApFtN2JyWdktbARtIpwB8BjwHbgY1pt43AnenxduCT6SrlWuDFiW6n\nmVUs0qKIOaVFOS2xM4A7muV/WAD8e0R8T9KPgO9I2gQ8BVyZ9r8buAwYAX4LfKrnUZvZQLQ9aJ9D\nUcDkNUm/AZ4YdByZlgCHBx1EhlrihHpirSVOmDrW34mIpd0eUNL30nFzHI6I9d2+10yUksR256zH\nXYJaYq0lTqgn1lrihLpinS2v7GpmVXMSM7OqlZLEpvw+ukLVEmstcUI9sdYSJ9QV66wUMSZmZtat\nUlpiZmZdcRIzs6oNPIlJWi/pibT+2JbpX9HXWG6UdEjSYx11Ra6bJmmlpPsk7ZX0uKRrS4xX0kJJ\nD0r6cYrzS6n+TEm7Upy3Sjox1Z+UtkfS86vaiLMj3iFJj0i6q/A4vcZfMtAkJmkI+AbNGmRnA1dL\nOnuAId0ETJ6gV+q6aUeAz0bEWcBa4Jr0uyst3peBdRFxDnAusD7djvYV4PoU5/PAprT/JuD5iHg3\ncH3ar03XAns7tkuNE7zGXyMiBlaADwA7OravA64bcEyrgMc6tp8AlqXHy4An0uN/Aq6ear8BxX0n\ncEnJ8QInAw8D76eZTb5g8v8DYAfwgfR4QdpPLcW3gubkXwfcRbMiS3Fxpvf8JbBkUl2x//b9LIPu\nTmatPTZgs1o3rQ2pK3MesIsC401dtD00K53sBH4BvBARR6aI5bU40/MvAovbiBP4GvA5YOIO5sWF\nxgl9WOOvVoP+opCstccKVUTskt4CfBf4dET8Ot2oP+WuU9S1Em9EjAHnSloE3AGcdZxYBhKnpI8A\nhyLiIUkXZcQy6H//nq/xV6tBt8RqWHvs6bReGqWtmybpBJoEdnNE3J6qi403Il4A7qcZw1skaeKP\naGcsr8WZnn878FwL4V0IfFTSL4Fv03Qpv1ZgnMAb1/ij+cPw2hp/Kaai/u37adBJ7EfA6nQF6ETg\nKpr1yEpS5LppappcW4G9EfHVUuOVtDS1wJD0ZuBDNAPn9wFXHCPOifivAO6NNJDTTxFxXUSsiIhV\nNP8P742Ij5UWJ3iNv6MMelCOZu2xn9OMk3xhwLHcAhwEXqX567WJZpzjHmBf+nla2lc0V1Z/AfwE\nWNNyrL9P0yV4FNiTymWlxQu8D3gkxfkY8Lep/l3AgzTrzv0HcFKqX5i2R9Lz7xrA/4OLgLtKjTPF\n9ONUHp84b0r7t2+r+LYjM6vaoLuTZmaz4iRmZlVzEjOzqjmJmVnVnMTMrGpOYmZWNScxM6va/wOT\nPpXlydZ2nwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'numpy.ndarray'>\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASQAAAD8CAYAAADe49kaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEgBJREFUeJzt3W+sZVV5x/Hvj5mhMDgj6qjRGVIg\nGlpCImMmVJ3EWMY2iAT6wiaQYNQ0GV9UxT+Ngb7o1HdtYgy+MKY3KJpAsThiaowiRKXEpAWHP2mB\ngZQiwoVRmBp0BJ2/v74459p7ZwbvPufutc/ad/8+yQn33Hvm2c/V4WGttfdaj2wTEVGDU2adQETE\nghSkiKhGClJEVCMFKSKqkYIUEdVIQYqIaqQgRUQ1UpAiohopSBFRjbUlgq6XfGaJwGP7uLBgdIAD\nheMDbCgc/3Dh+C8Wjr+mcHyAo4XjHywY+wXsl7SSCG+S/FLDz+6D79m+ZCXXa6JIQToT2Fki8Nhn\n+LeC0QHuKhwf4F2F488Xjn9v4fivLhwf4BeF4z9ZMPbciiO8BHy44Wf/Hjat+IINFClIEVE/UV8B\nqC2fiOjIKcDps07iOClIEQMlYN2skzhOClLEQGXKFhHVyAgpIqqREVJEVKPGEVKjJ7UlXSLpMUmP\nS7q2dFIRUd7CXbYmr+VI+rKk5yQ9tOh7r5Z0p6T/Hv/zVU1yWu5Ca4AvAO8BzgeuknR+gxwjomIL\nI6Qmrwa+Ahz/JPe1wPdtvxn4/vj979VkhHQR8LjtJ2wfAr4GXNEsx4io2dqGr+XYvpsTH32/Avjq\n+OuvAn/RJJ/lbAaeXvR+HviT4z8kaSfjHSOvbBA0ImargzWk19veB2B7n6TXLfcHmhSkk23gO6F3\nku05xhts3iilt1JE5Sa8y7ZJ0p5F7+fG/863qkk+88BZi95vAZ5tO5GI6NaEW0f229424SV+LukN\n49HRG4DnmuS0nB8Db5Z0jqRTgSuBb02YWERUpuVF7ZP5FvCB8dcfAP51uT+w7AjJ9hFJHwG+x+iQ\nmi/bfnj6HCOiBm0+GCnpFkZn6mySNA/sAv4BuFXSXwFPAX+5XJxG+dj+DvCdqbONiOq0uaht+6qX\n+dGOSeLkSe2IgcrWkYioRo1bR1KQIgZK5IC2iKiEgHVNK8CRkpn8vxSkiIGSYG0KUkTUQIJ1XXSb\nmkCRgrSPC4u2KtpVeLfcZ/hl0fgjpdsUbSkcv7TSbZagm1ZLpfzBiiNMNELqSGXpRERXJFi38rrW\nqhSkiKGq8EGkytKJiM6kIEVEVSqrAJWlExGdEaPt8hVJQYoYqkzZIqIaoo2nB1qVghQxVBWOkJq0\nQTqh31JErAILBamNtiMtaXKE7Vc4sd9SRKwGaxq+OtLkCNu7JZ1dPpWI6FSFU7bK0omIzqzmgrS4\nUeTSrkkRUaXVfJdtcaNIaWsaRUbUbjWPkCKiZyosSE1u+98C/DtwnqT5cY+liOi7ha0jPbvL9nL9\nliKizyocIVWWTkR0ZjUvakdEz2SEFBHVSEGKiKpUVgEqSyciOpMD2iKiGsOZsh0A7ioTmvJ900r3\nfQP4DDcWv0ZZpfu+3Vo4fhdK9n1rYWiTu2wRUY3hjJAionoVFqQmB7RFxGrU4omRkj4h6WFJD0m6\nRdJp06SUghQxZC3sZZO0GfgYsM32BeM/ceU06VQ2YIuIzrQ7ZVsLnC7pMLAeeHbaIBExRJPdZdsk\nac+i93PjM9Cw/YykzwJPAb8B7rB9xzQppSBFDNVkI6T9tredNIz0KuAK4BzgBeDrkq62fdOkKWUN\nKWKo2lvUfjfwE9vP2z4M3Aa8Y5qUmhzQdpakH0raO15Fv2aaC0VEZdo7oO0p4G2S1ksSsAPYO01K\nTQZsR4BP2b5f0gbgPkl32n5kmgtGRCVaWtS2fY+k3cD9jOrFA4zP159UkxMj9wH7xl8fkLQX2Ayk\nIEX0mYCpnhY6ke1dwK6VxpmoPo4bRm4F7lnphSNixvq821/SK4BvAB+3/auT/HxRX7bXtpReRBRT\n4daRRulIWseoGN1s+7aTfWZpX7Y3pS9bRB/0rSCNV82/BOy1/bnyKUVEJyqcsjV5Dmk78H7gYkkP\njl+XFs4rIkprcXNtW5rcZfsRo9QjYjXJAW0RUY2+LmpHxCqUghQR1UhBioiqVHaXLQUpYqgyQoqI\nauQuW0RUYzgjpA3Au8qEBmC+YOxumjju4kNF45duplmyEejIJwvHh9J/j+DegrGPrjzEcApSRFQv\nBSkiauLcZYuIGvgUONTSAW1tSUGKGCgLjqxp2ufjWNFcFqQgRQyUJY6ubVoCDhXNZUEKUsSAHV1T\n1yJSClLEQBlxtLK9I01OjDwNuJvRM51rgd3jDgMR0WNGHOlbQQIOAhfb/vX4bO0fSfqu7f8onFtE\nFGTEocr2jjQ5MdLAr8dv141fOcQ/oud6OWUDkLQGuA94E/AF2+nLFrEK1FaQGj2EYPuo7QuBLcBF\nki44/jOSdkraI2kP/G/beUZEyxbWkJq8utL0qSgAbL/AaFflJSf52Zztbba3wWtaSi8iShlN2dY2\nenWlyV221wKHbb8g6XTg3cA/Fs8sIooaLWqfOus0lmhS+t4AfHW8jnQKcKvtb5dNKyJKM/Tvtr/t\n/wS2dpBLRHRKnU7Hmqgrm4joTG9v+0fE6lRbQZroLltErB4LI6Qmr+VIOlPSbkmPStor6e3T5JQR\nUsRAGXGwva0jnwdut/0+SacC66cJkoIUMVBtrSFJ2gi8E/gggO1DTHmAUqZsEQPV4pTtXOB54EZJ\nD0i6QdIZ0+SUghQxYBNsHdm0sDVs/Nq5KMxa4K3AF21vBV4Erp0mn0JTtsOU7Xm1pWDsbpTum7aL\nVxaNX77vW+meadDvv0dTDUCW8GTPIe0fbQs7qXlgftGm+93UVZAionZtrSHZ/pmkpyWdZ/sxYAfw\nyDSxUpAiBmp0l621vWwfBW4e32F7AqZrzZyCFDFQE07Zfn8s+0Hg5aZ0jaUgRQxYbU9qpyBFDFT2\nskVENVKQIqIaLW8daUXjgjQ+oG0P8Izty8qlFBFd6PsI6RpgL7CxUC4R0bHaClKjrSOStgDvBW4o\nm05EdKXGriNNR0jXA58GNhTMJSI61OZzSG1ZdoQk6TLgOdv3LfO5RX3ZftFaghFRTlsHtLWlSXnc\nDlwu6VLgNGCjpJtsX734Q7bngDkA6YK02o6oXI1tkJYdIdm+zvYW22cDVwI/OL4YRUT/9HkNKSJW\nmRrXkCbKxvZdjFppR8QqUNtt/7rKY0R0pu8PRkbEKrKwhlSTFKSIgRrdZevpXraIWF0yZYuIqqQg\nRUQVsoYUEdXo/XNIzb0I3FsmdCe66Nd1V9Ho/e/7dmPR+N0o+fdo3Yoj1Lh1pK7yGBGdyZQtIqoy\nkClbRNQut/0johopSBFRlawhRUQVjnFKto5ERD0yZYuIKvR2DUnSk8AB4ChwxPa2kklFRHmm32tI\nf2p7f7FMIqJjg9k6EhG16+2UjdHo7g5JBv5p3PJoCUk7gZ2jd69pK7+IKMSIgy3uZZO0BtgDPGP7\nsmliNC1I220/K+l1wJ2SHrV99+IPLO3Ldk76skVUrsBu/2uAvcDGaQMs25cNwPaz438+B3wTuGja\nC0ZEPdrqXCtpC/Be4IaV5NOklfYZkjYsfA38OfDQSi4aEbO3sIbUUivt64FPA8dWklOT8drrgW9K\nWvj8P9u+fSUXjYjZM+LoscaL2psk7Vn0fm5hLVnSZcBztu+T9K6V5LRsQbL9BPCWlVwkIurjY+Lg\nbxtvHdn/e54/3A5cLulS4DRgo6SbbF89aU657R8xULY4emTlt/1tXwdcBzAeIf3NNMUIUpAihsu0\nUpDalIIUMVC2OHK43YJk+y5WcGB8ClLEYIljR+sqAXVlExHdMZApW0RU4Zjgt3WVgELZrAFeXSY0\nUL7n262F4wN8snD8+aLRS/dN28WHisaH8r3ryvbeO9BOmCPthGlLXeUxIrozOhCpKilIEUOVghQR\n1TBweNZJLJWCFDFUBg7OOomlUpAihipTtoioRgpSRFQjBSkiqlFhQWp0hK2kMyXtlvSopL2S3l46\nsYjowJGGr440HSF9Hrjd9vsknQqsL5hTRHThGPDbWSex1LIFSdJG4J3ABwFsHwIOlU0rIorr6ZTt\nXOB54EZJD0i6YXzY/xKSdkraMzp391etJxoRLVsoSBVN2ZoUpLXAW4Ev2t4KvAhce/yHbM/Z3jY6\nd3fqtkwR0ZWeFqR5YN72PeP3uxkVqIjou74VJNs/A56WdN74WzuAR4pmFRHlVThCanqX7aPAzeM7\nbE9AB4fVRERZx4DfzDqJpRoVJNsPAi/Xkyki+sjA0VknsVSe1I4Ysspu+6cgRQxVhc8hpSBFDFUK\nUkRUo49bRyJiFcsIKSKqkClbRFRjOIf8HwV+USY0ULYJZVfKNnKELYXjl1W+iSPs4pVF45f9HTas\nPESeQ4qIamTKFhHVMP3cOhIRq1CmbBFRjQqnbI0O+Y+IVail40cknSXph+MGIA9LumbalDJCihiq\n9m77HwE+Zft+SRuA+yTdaXvic9NSkCKGrIU1JNv7gH3jrw9I2gtsZoqDHJt0HTkP+JdF3zoX+Dvb\n1096sYioyGR72TaNGnj8zpztueM/JOlsYCtwz/E/a2LZgmT7MeDC8cXWAM8A35zmYhFRkcmmbPtH\nDTxenqRXAN8APm57qtZDk07ZdgD/Y/un01wsIirS4m1/SesYFaObbd82bZxJC9KVwC3TXiwiKtPC\nbX9JAr4E7LX9uZXEanzbf3zA/+XA11/m54saRR5YSU4R0YX2uo5sB94PXCzpwfHr0mlSmmSE9B7g\nfts/P9kPxwtccwDSOZ4mmYjoUEsHtNn+EaCVR5qsIF1FpmsRq0eFT2o3KkiS1gN/Bny4bDoR0ak+\nFiTbLwGvKZxLRHRpOAe0RUT1sts/IqrR1zWkiFiFjpED2iKiIpmyRUQ1KntiMAe0RUQ1UpAiohqF\npmwHgSfLhO5EF33f7u3gGiWV7vt2V+H45Xu/lez7dsJBRKtE1pAiBqu+22wpSBGDVd+j2ilIEYNV\n35ORKUgRg5URUkRUIwUpIqphsqgdEZWobw2p0YORkj4xbpH7kKRbJJ1WOrGIKG1hytbk1Y1lC5Kk\nzcDHgG22LwDWMOo+EhG91t4p/21pOmVbC5wu6TCwHni2XEoR0Y36FrWXHSHZfgb4LPAUo/7dv7R9\nR+nEIqK0+kZITaZsrwKuAM4B3gicIenqk3xuUV+2l9rPNCJatrB1pMmrG00Wtd8N/MT287YPA7cB\n7zj+Q7bnbG8b9f9e33aeEdG6+ha1m6whPQW8bdwK6TfADmBP0awioiN13fZftiDZvkfSbuB+Rtk/\nwOo9/SBiQOpb1G7al20XsKtwLhHRqZ4WpIhYjep7UjsFKWKwckBbRFQjU7aIqEambBFRjYyQIqIa\nGSFFRDXqW9SW3X4vXUnPAz+d4I9sAva3nkh3kv/s9f13mDT/P7T92pVcUNLt4+s2sd/2JSu5XhNF\nCtLESUh7Rnvg+in5z17ff4e+59+WtNKOiGqkIEVENWopSH3frJv8Z6/vv0Pf829FFWtIERFQzwgp\nImK2BUnSJZIek/S4pGtnmcs0JJ0l6YeS9o7bRF0z65ymIWmNpAckfXvWuUxK0pmSdkt6dPz/w9tn\nndMk0mJsqZkVJElrgC8A7wHOB66SdP6s8pnSEeBTtv8YeBvw1z38HQCuAfbOOokpfR643fYfAW+h\nR79HWoydaJYjpIuAx20/YfsQ8DVGzQR6w/Y+2/ePvz7A6F+GzbPNajKStgDvBW6YdS6TkrQReCfw\nJQDbh2y/MNusJrbQYmwtaTE204K0GXh60ft5evYv82KSzga2AvfMNpOJXQ98mtE+gr45F3geuHE8\n5bxB0hmzTqqptBg70SwLkk7yvV7e8pP0CuAbwMdt/2rW+TQl6TLgOdv3zTqXKa0F3gp80fZW4EWg\nN2uRTVuMDcksC9I8cNai91vo4XBV0jpGxehm27fNOp8JbQcul/QkoynzxZJumm1KE5kH5m0vjEp3\nMypQfdGoxdiQzLIg/Rh4s6RzJJ3KaDHvWzPMZ2KSxGj9Yq/tz806n0nZvs72FttnM/rf/we2e/Nf\naNs/A56WdN74WzuAR2aY0qR+12Js/HdpBz1alC9hZseP2D4i6SPA9xjdXfiy7Ydnlc+UtgPvB/5L\n0oPj7/2t7e/MMKeh+Shw8/g/ak8AH5pxPo2lxdiJ8qR2RFQjT2pHRDVSkCKiGilIEVGNFKSIqEYK\nUkRUIwUpIqqRghQR1UhBiohq/B/Ynr6gQJpfUwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "## lattice and material parameters\n",
    "a = 1;\n",
    "radius = 0.2*a; #0.4;\n",
    "e_r = 12;\n",
    "\n",
    "#generate irreducible BZ sample\n",
    "T1 = 2*np.pi/a;\n",
    "T2 = 2*np.pi/a;\n",
    "\n",
    "## Specify number of fourier orders to use:\n",
    "N = 1; M = 1;\n",
    "NM = (2*N+1)*(2*M+1);\n",
    "\n",
    "# ============== build high resolution circle ==================\n",
    "Nx = 512; Ny = 512;\n",
    "A = e_r*np.ones((Nx,Ny));\n",
    "ci = int(Nx/2); cj= int(Ny/2);\n",
    "cr = (radius/a)*Nx;\n",
    "I,J=np.meshgrid(np.arange(A.shape[0]),np.arange(A.shape[1]));\n",
    "\n",
    "dist = np.sqrt((I-ci)**2 + (J-cj)**2);\n",
    "A[np.where(dist<cr)] = 1;\n",
    "\n",
    "#visualize structure\n",
    "plt.imshow(A);\n",
    "plt.colorbar();\n",
    "plt.show()\n",
    "\n",
    "## =============== Convolution Matrices ==============\n",
    "\n",
    "## fft alone\n",
    "Af = (1 / np.prod(N)) * np.fft.fftshift(np.fft.fft2(A));\n",
    "\n",
    "\n",
    "E_r = cm.convmat2D(A, N,M)\n",
    "print(type(E_r))\n",
    "plt.figure();\n",
    "plt.imshow(abs(E_r), cmap = 'jet');\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ -11.        +0.00000000e+00j    2.78795942-1.13686838e-13j\n",
      "    54.49429876+0.00000000e+00j ... -179.77229926+1.43884904e-13j\n",
      "    54.49429876+0.00000000e+00j    2.78795942+1.13686838e-13j]\n",
      " [   2.78795942-3.77054237e-16j  -19.87525009+6.25855302e-14j\n",
      "    19.54104782+5.39070621e-14j ...   62.54158308-6.48181674e-14j\n",
      "    19.54104782-6.24097146e-14j  -19.87525009-2.62828830e-14j]\n",
      " [  54.49429876-3.82805798e-15j   19.54104782+4.55334635e-14j\n",
      "  -150.11657726-7.96018224e-14j ...  174.91110779-2.63467932e-14j\n",
      "  -150.11657726+1.25818226e-13j   19.54104782-3.91610824e-14j]\n",
      " ...\n",
      " [-179.77229926-3.57932703e-14j   62.54158308+9.77834630e-14j\n",
      "   174.91110779+5.46843590e-14j ... -306.47861514-1.16587254e-13j\n",
      "   174.91110779-5.05109893e-15j   62.54158308-6.59190484e-14j]\n",
      " [  54.49429876+1.71368991e-14j   19.54104782+3.91610824e-14j\n",
      "  -150.11657726-1.25818226e-13j ...  174.91110779-1.46240502e-15j\n",
      "  -150.11657726+7.96018224e-14j   19.54104782-4.55334635e-14j]\n",
      " [   2.78795942-2.87799670e-15j  -19.87525009+2.62828830e-14j\n",
      "    19.54104782+6.24097146e-14j ...   62.54158308-3.24118763e-14j\n",
      "    19.54104782-5.39070621e-14j  -19.87525009-6.25855302e-14j]]\n"
     ]
    }
   ],
   "source": [
    "np.set_printoptions(precision = 8)\n",
    "print(Af)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Matlab answer AT LEAST agrees in terms of the magnitude (abs) but the real and imaginary parts do appear to have disagreement\n",
    "ans =\n",
    "\n",
    "   1.0e+02 *\n",
    "\n",
    "  Column 1\n",
    "\n",
    " -0.110000000000000 + 0.000000000000000i  <br>\n",
    "  0.027877494914885 + 0.000342125507580i <br>\n",
    "  0.544778860922760 + 0.013373570389723i <br>\n",
    "\n",
    "  Column 2\n",
    "\n",
    "  0.027877494914876 + 0.000342125507579i <br>\n",
    " -0.198692640333564 - 0.004877630543374i <br>\n",
    "  0.195278065014253 + 0.007192517034929i<br>\n",
    "\n",
    "  Column 3\n",
    "\n",
    "  0.544778860922757 + 0.013373570389724i <br>\n",
    "  0.195278065014249 + 0.007192517034928i <br>\n",
    " -1.499357552702474 - 0.073658713242330i <br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
